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Structural, elastic and thermal properties of cementite (FeaC) were studied using a Modified 
Embedded Atom Method (MEAM) potential for iron-carbon (Fe-C) alloys. Previously developed 
Fe and C single element potentials were used to develop a Fe-C alloy MEAM potential, using a 
statistically-based optimization scheme to reproduce the heat of formation of Fe-C alloys in L12, 
Bi, and cementite structures, as well as the interstitial energies of C in bcc Fe. The stability of ce- 
mentite was investigated by molecular dynamics simulations at high temperatures. The nine single 
crystal elastic constants for cementite were obtained by computing total energies for strained cells. 
Polycrystalline elastic moduli for cementite were calculated from the single crystal elastic constants 
of cementite. The formation energies of (001), (010), and (100) surfaces of cementite were also 
calculated. The melting temperature and the variation of specific heat and volume with respect to 
temperature were investigated by performing a two-phase (solid/liquid) molecular dynamics sim- 
ulation of cementite. The predictions of the potential are in good agreement with first-principles 
calculations and experiments. 

PACS numbers: 61.50.Lt, 62.20.de, 61.72.jj, 68.35.Md, 71.15.Pd 



I. INTRODUCTION 



Steel alloys are the most widely used structural materi- 
als due to their abundance, all-purpose applicability and 
low cost. The main carbide in steel alloys is cementite, 
which forms as a precipitate. Cementite has a direct im- 
pact on the mechanical, structural, and thermal proper- 
ties of steel. Therefore, the ability to describe and predict 
properties of cementite at the nanoscale is essential in 
the study and design of new steels. Atomistic simulation 
methods, such as molecular dynamics or Monte Carlo 
simulations, offer an efficient and reliable route to inves- 
tigating nanoscale mechanics pertaining to cementite in 
steel alloys. Each of these methods requires accurate in- 
teratomic potentials to find the energy of the system un- 
der investigation. However, first-principles calculations- 
albeit the most reliable-are incapable of simulating the 
number of atoms required for realistic calculations due to 
unreasonable memory and processing-time requirements. 



Therefore, semi-empirical potential methods are being 
explored as a suitable alternative. 

Among the spectrum of semi-empirical formulations, 
the Modified Embedded Atom Method (MEAM), orig- 
inally proposed by Baskes et has been shown to 
accurately predict properties of most crystal structures, 
such as bcc, fee, hep, and even diatomic gases, in good 
agreement with experiments or first-principles calcula- 
tions. MEAM is extended from the Embedded Atom 
Method (EAMJPl to include the directionality of bonds. 
In the original MEAM formalism, only the first-nearest 
neighbor (INN) interactions were considered.'^' Lee and 
BaskePI later extended the original formalism to include 
the screened second-nearest neighbor (2NN) interactions. 
Further details of the MEAM formalism can be found in 
Ref.lIlandU 

One of the commonly used 2NN MEAM potential for 
the Fe-C system developed by Byeong-Joo Lec^ is de- 
signed to predict the interactions of interstitial C atoms 
with defects such as vacancies. According to Fang et 
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al.J^ Lee's potential predicts a negative heat of forma- 
tion for cementite. Lee's potential also predicts that ce- 
mentite is only stable up to a temperature of 750 KP 
Experimentally, however, cementite is metastable with a 
positive heat of formationP Moreover, according to the 
Fe-C phase diagram,'^ cementite is stable up to 1420 K. 
Among recent interatomic potential^SHUfor the Fe-C sys- 
tem, EAM potentials by Lau et alPand Ruda et al.^^ and 
the analytical bond order potential (ABOP) by Henriks- 
son et all^^l all promise to predict properties of cementite 
reasonably well. In the potentials by Lau et alW and 
Ruda et al.,1^ however, the single element potential for 
C does not predict properties of both graphite and di- 
amond well. This is due to the limited ability of EAM 
to describe the bare C-C interaction successfully.^^ We 
note that a successful interatomic potential for an alloy 
system should not only predict the properties of the al- 
loy correctly, but it should also predict the properties of 
the individual alloying elements in their natural crystal 
structures accurately. The ABOP by Henriksson et al.l^ 
predicts properties of cementite as well as Fe and C accu- 
rately; however, ABOPs are not applicable to simulations 
involving interfaces and surfaces; Also AB OPs a re re- 
stricted to considering only INN interactions !^^ * ^^ ! Some 
of the more recent potentials for the Fe-C system are im- 
plemented using in-house developed molecular dynamics 
codes, which limits the potentials' transferability. 

In the present work, we develop a 2NN MEAM po- 
tential for the Fe-C alloy system that accurately predicts 
the structure and properties of cementite. Our Fe-C alloy 
potential is based on previously developed 2NN MEAM 
potentials for F^^^ and CC^I in their pure forms. The C 
MEAM potential predicts both diamond and graphite as 
minimum energy structures with almost degenerate en- 
ergies. Using the Fe and C single element potentials, the 
alloy potential for Fe-C was constructed by fitting to the 
heat of formation of Fe-C alloys in the Bi and L12 struc- 
tures, and cementite, as well as the interstitial energies 
of C in bcc Fe. 



II. METHODS 

For all atomistic simulations described in this p aper, 
we have used MEAM as implemented in LAMMPS,^^ 
the classical molecular dynamics simulation code from 
Sandia National Laboratories. Some of the reference 
data required for potential construction and validation 
are not readily available from experiments. In these 
cases, we performed the first-principle s calc ulations us- 
ing Density Functional Theory (DFT)2212U and Projec- 
tor Augmented Wave (PAW) pseudopotentials.'^ Elec- 
tron exchange and correlation were treated with Gener- 
alized Gradient Ap proximation (GGA) as parameterized 
by Perdew et al.'^^lBrillouin zone sampling was performed 
using the Monkhorst-Pack scheme,'2Sl with a Fermi-level 
smearingof 0.2 eV applied using the Methfessel-Paxton 
method.ES Geometric optimizations were carried out us- 



ing the conjugate gradient minimization method.'^fll The 
constructed Fe-C alloy potential was optimized through 
a statistically-based optimization scheme to reproduce 
heat of formation of Fe-C alloys in the Bi and L12 struc- 
tures, and cementite, as well as the interstitial energies 
of C in bcc Fe. 

The MEAM potential was constructed by adjusting the 
potential parameters to reproduce several selected target 
values using a statistically-based optimization scheme de- 
veloped by Tschopp et This scheme uses Latin Hy- 
percube Sampling (LIIS),^l-a stratified random sampling 
approach-to sample the A'^-parameter space, where N is 
the number of potential parameters sensitive to the tar- 
get values chosen. A sensitivity analysis was completed 
at the beginning of the optimization to determine param- 
eters of the potential model most sensitive to the targets. 
Next, random sampling of the parameter space was done 
using LHS, with 4000 different potential parameter com- 
binations. Then analytical models were generated to pro- 
vide a representation of the correlation between poten- 
tial parameters and target values. Once the models were 
generated, a multi-objective optimization was executed 
to determine the best possible values for the parameters. 

III. MEAM POTENTIALS 

A. Single element potentials 

The single element MEAM potential parameters are 
presented in Table |l] The parameters for Fe are from 
the MEAM potential developed by Lee et al.,1^ and the 
parameters for C are from Uddin et al.^^ 

1. Energy vs. volume curves 

Energy variation with respect to volume or nearest 
neighbor distance is considered an important test of va- 
lidity of interatomic potentials. Here we present the en- 
ergy vs. volume curves generated by the single element 
potential for Fe and energy vs. nearest neighbor dis- 
tance curves generated by the single element potential 
for C. Fig. [1] shows the energy vs. volume curve for bcc 
Fe compared to the curves generated by DFT calcula- 
tions and experimental data. It is well known that DFT 
overestimates the cohesive energy.'^H! Therefore, the DFT 
curve is shifted vertically by a constant amount to the 
experimental cohesive energy at the equilibrium volume 
to aid the comparison of the curves. Due to overbind- 
ing, the DFT's prediction for the equilibrium volume is 
underestimated.'22l Therefore, the DFT curve sits to the 
left of the experimental curve. The experimental curve 
was constructed from Murnaghan's equation of stat^SSEI] 
using the experimental bulk modulus, cohesive energy, 
and atomic volume listed in Table |lll We also tested the 
stability of Fe in several different crystal structures in- 
cluding bcc, fee and hep structures as shown in Fig. [2] 
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TABLE I. Set of the MEAM potential parameters for pure Fe and C. The parameters include: cohesive energy per atom Ec 
(eV), equilibrium nearest neighbor distance of the reference structure re (A), cutoff radius rcut (A), density scaling factor po, 
embedding energy scaling factor A, exponential decay factor for the universal energy a, additional cubic term in the universal 
energy equation as, atomic density exponential decay factors I3^'^~^\ weighting factors for atomic densities t'^"^', and angular 
screening parameters Cmax and Cmin- The bcc and diamond lattices are chosen as the reference structures for Fe and C, 
respectively. 
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FIG. 1. Comparison of energy vs. volume curves for Fe in 
its most stable bcc structure. The solid curve is constructed 
from experimental values in Table |ll] For comparison, the 
DFT curve is shifted vertically to match the experimental 
cohesive energy of Fe at the equilibrium volume. 
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FIG. 2. Comparison of energy vs. volume curves for Fe 
in bcc, fee and hep crystal structures. For ease of compari- 
son, the DFT curves are shifted vertically by a constant value 
equal to the difference between experimental and DFT cohe- 
sive energies of Fe in bcc at equilibrium volumes. 



The Fe MEAM potential correctly predicts that bcc is 
the most stable structure, as observed in experiment and 
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FIG. 3. Comparison of energy vs. nearest neighbor distance 
curves for C in diamond and graphite structures. The solid 
curve is constructed from experimental values in Table [III F^i" 
comparison, the DFT curve is shifted vertically to the exper- 
imental cohesive energy at the equilibrium nearest neighbor 
distance. 




FIG. 4. Cohesive energy of graphite as a function of the c/a 
ratio. Energy at zero is set to the minimum energy predicted 
by MEAM. 



the first-principles methods. MEAM predicts that fee 
and hep Fe are much closer in energy and have a larger 
volume than that calculated from DFT. 
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C single element MEAM potential predicts both di- 
amond and graphite to be stable structures. Energy 
vs. nearest neighbor distance curves for diamond and 
graphite are shown in Fig. [3j The experimental curves 
were constructed from the Rose's equation of state^^^' us- 
ing the experimental bulk modulus, cohesive energy and 
nearest neighbor distance, as listed in Table |TTj MEAM 
predictions for diamond and graphite are in good agree- 
ment with experiment. MEAM predicts almost degen- 
erate cohesive energies for graphite and diamond, while 
DFT predicts graphite to be ~ 0.1 eV more stable than 
diamond. For graphite, DFT predicts a first-nearest 
neighbor (INN) distance in good agreement with ex- 
periment, while MEAM predicts a INN distance 3% 
greater than the experimental value. MEAM optimized 
the c/a ratio of the graphite structure to 3.35. This does 
not agree well with the experimental value of 2.725:-^ 
This disagreement is due to the incorrect prediction of in- 
terlayer interactions of graphite, which are dominated by 
van Der Waals forces, that are not described by MEAM. 
However, the dependence of cohesive energy on the c/a 
ratio is small. As illustrated in Fig. |4j the difference 
in cohesive energy of graphite between the experimental 
and MEAM c/a ratio is ~ 4 meV/atom. In construct- 
ing the energy vs. nearest neighbor distance curves for 
graphite, the inter-planar distance was scaled with the 
lattice constant. The experimental ratio was used in the 
generation of the DFT curve, while the MEAM curve was 
constructed with the predicted c/a ratio. 



2. Single element material properties 

The cohesive energy, equilibrium lattice constants and 
bulk moduli for bcc Fe, graphite, and diamond were de- 
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termined by fitting Murnaghan's equation of state 
the energy vs. volume curves generated by MEAM. The 
single element material properties compared to experi- 
mental values are given in Table lllj 



TABLE II. Material properties predicted by the single 
element MEAM potentials. Ec is the cohesive energy 
(eV/atom); a and c are the equilibrium lattice constants (A); 
B is the bulk modulus (GPa); and fio is the equilibrium 
atomic volume (A^/atom). Experimental data are given in 
parentheses. Experimental values for equilibrium atomic vol- 
ume were calculated from the experimental lattice parame- 
ter(s). 



Property 



B 
a 

c 
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177 (16^''^ 
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^ Philipsen et al.I^Hl 
^ Wang et alPl 

" Ref.|35|and|36]as reported by Fahy et all^ll 
Ref.|30]and|38]as reported by Yin et al.Ell 



L12, and cementite structures. The final parameters of 



the alloy potential are presented in Table III 



Table|IV]lists the MEAM prediction of the target prop- 
erties used for fitting of the potential compared to DFT 
and experimental values, and to the interatomic poten- 
tials by Ruda et al.^^ and Henriksson et al.'^ The or- 
ders of stability of Bi, L12, and cementite structures are 
predicted correctly compared with DFT results. MEAM 
predicts the heat of formation of Bi and L12 structures in 
good agreement with DFT values. MEAM prediction of 
the heat of formation of cementite is also in good agree- 
ment with DFT and experimental data. The MEAM- 
predicted interstitial energies of C in tetrahedral and oc- 
tahedral positions of bcc Fe have the same order of mag- 
nitude as in DFT predictions. Also like DFT, MEAM 
predicts the octahedral interstitial to be more stable than 
the tetrahedral interstitial. 



B. MEAM potential for the Fe-C alloy 

The MEAM potential parameters for the Fe-C alloy 
potential were constructed by fitting the heat of forma- 
tion of Bi, L12, and cementite, and to the interstitial en- 
ergies of C in bcc Fe. Fe3C in the L12 structure and FeC 
in the Bi structure are hypothetical crystal structures 
that do not exist in nature. Therefore, properties of Fe-C 
in L12 and Bi crystal structures were obtained by first- 
principles calculations. We also used the first-principles 
calculations to obtain the interstitial energies of C in the 
bcc Fe lattice at octahedral and tetrahedral positions. 
The potential was constructed using a statistically-based 
optimization scheme developed by Tschopp et al.'^Sl More 
weight was given to the set of target values that repro- 
duces the correct order of stability of the two C intersti- 
tial energies in bcc Fe and the heat of formation of Bi, 



1. Energy vs. volume curves for B\ and L12 structures 

The cohesive energy of Fe-C in the Bi and L12 crystal 
structures as a function of the atomic volume is shown 
in Figs. [5] and |6] respectively. For the Bi structure, 
MEAM predicts an atomic volume ~ 11% less, and a 
bulk modulus 0.3% less than the DFT results. The 
MEAM prediction for the L12 structure gives an atomic 
volume ~ 11% greater, and a bulk modulus 35% less than 
DFT. As mentioned earlier, DFT overestimates the co- 
hesive energy. Therefore, to aid the comparison in these 
figures, the DFT curve is shifted vertically by a constant 
amount to the MEAM-predicted cohesive energy at the 
equilibrium nearest neighbor distance. 
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TABLE III. The MEAM potential parameters for the Fe-C alloy system. The parameters include: heat of formation of 
the reference structure A [eV] , equihbrium nearest neighbor distance of the reference structure re [A] , cutoff radius rcut [A] , 
embedding energy scaling factor A, exponential decay factor for the universal energy a, additional cubic term in the universal 
energy equation as, and angular screening parameters Cmax and Cmin. The triplet {a,b,c) represents the configuration with c 
atom in between a and b atoms. The Bi lattice is chosen as the reference structure. 
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TABLE IV. Materials properties of Fe-C alloys used for the MEAM potential construction. Heat of formation per atom 
(eV/atom) of Fe-C alloys in Bi, L12, and cementite structures; the interstitial energies of C in bcc FefeV). The predictions 
from the present work (MEAM) are compared with those from DFT and other interatomic potentialj^^^^. Experimental data 
are given in parentheses. 
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FIG. 5. Comparison of the energy vs. volume curves of Fe-C 
alloy system in the Bi structure. DFT curve is shifted ver- 
tically to the MEAM-predicted cohesive energy at the equi- 
librium nearest neighbor distance to aid the comparison with 
the MEAM curve. 
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FIG. 6. Comparison of the energy vs. volume curves of Fe-C 
alloy system in the L12 structure. DFT curve is shifted verti- 
cally to the MEAM-predicted cohesive energy at the equilib- 
rium nearest neighbor distance aid the comparison with the 
MEAM curve. 



2. Elastic constants of FeC in the Bi crystal structure 



Elastic constants of Fe-C in the Bi crystal structure 
were calculated using the Fe-C MEAM potential and 
listed in Table |V|in comparison with DFT calculations. 
They were calculated using the deformation matrix pre- 
sented in Jiang et al.^^ In linear elastic theory, defor- 
mation energy is a function of strain. Distortion ener- 
gies (AE) calculated for strains {6) equal to ±2% and 
±0.1% were fitted to AE = fc2(5^ + hS^. DFT calcu- 
lations were performed for S = ±2% for comparison. 



The single-crystal elastic constants were obtained using 
the relationships for the quadratic coefficient fc2 hsted in 
Jiang et al. The present MEAM's result for cn compares 
reasonably well with the DFT result. C12 is predicted 
at a lower value than DFT, but it is in the same order 
of magnitude. MEAM prediction of C44 is significantly 
larger than the DFT result. 
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TABLE V. Elastic moduli (GPa) of FeC in the Bi structure. 
The predictions from the present work (MEAM) are compared 
with those from DPT. 



Property 


DFT 




MEAM 




3 = ±2% 


S = ±2% 


5 


= ±0.1% 


Cll 


601 


570 




566 


Cl2 


589 


213 




213 


C44 


83 


391 




389 



IV. STRUCTURAL AND ELASTIC 
PROPERTIES OF CEMENTITE 

We tested the present MEAM potential further by 
computing the structural properties of cementite includ- 
ing the equilibrium lattice parameters, the equilibrium 
volume per atom, and the heat of formation. These 
properties are presented in Table |VI| with comparison to 
DFT^xperiment, and interatomic potentials by Ruda 
et al!^ and Henriksson et As noted in Sec. IlIIB 



MEAM prediction of the heat of formation of cemen- 
tite is in excellent agreement with DFT and experimen- 
tal data. Henriksson's potential also predicts a value in 
good agreement with DFT and experiment, while Ruda's 
potential predicts a much larger value. Lattic e cons tants 
of the present MEAM and literature potentialiSMl agree 
well with experiment, while DFT predicts lower values. 
As a test of validity, the variation of cohesive energy 
with volume was calculated. Fig. [7] compares the en- 
ergy vs. volume curves for cementite generated by the 
present MEAM potential with DFT and experimental 
curves. During volume variation of cementite, the ra- 
tios between a, b and c lattice parameters were held con- 
stant. As noted before, DFT overestimates the cohesive 
energy and underestimates the volume. Therefore, the 
DFT curve sits to the left of the experimental curve, and 
it is shifted vertically to the experimental cohesive energy 
at equilibrium volume to aid the comparison with the ex- 
perimental curve. The experiment al cur ve was generated 
by Murnaghan's equation of stat^^SEU ^jth the exper- 
imental bulk modulus, volume, and cohesive energy.l^S! 
The experimental single-crystal bulk modulus of cemen- 
tite has not yet been determined; therefore, the polycrys- 
talline bulk modulus of cementite was used to generate 
the experimental curve. 



A. Single-crystal elastic properties 

The elastic moduli of cementite were calculated and 
compared to DFT and the interatomic potentials by 
Ruda et ali^ and Henriksson et al.^^ as presented in Ta- 
ble Ivnl We used both S = 0.5% and 5 = 0.1% that 
produced the same result. These results show that the 
present MEAM potential for Fe-C alloy predicts cemen- 
tite to be stable (positive elastic constants) and their val- 
ues are reasonably close to those predicted by DFT. How- 
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FIG. 7. Comparison of energy vs. volume curves for ce- 
mentite. The sohd curve is constructed from experimental 
values of the cohesive energy, polycrystalline bulk modulus, 
and equilibrium volume of cementite. For comparison, the 
DFT curve is shifted vertically to the experimental cohesive 
energy at the equilibrium volume. 



TABLE VI. Material properties of cementite: the heat 
of formation AHf (eV/atom); the equilibrium volume Qq 
(A"^/atom); the lattice parameters a, b, and c (A). The pre- 
dictions from the present work (MEAM) are com pared with 
those from DFT and other interatomic potentialS^^. Ex- 
perimental values are given in parentheses. 

Property DFT (Expt.) MEAM Ruda Henriksson 



no 
a 
b 
c 



0.01 (0.0 
8.49 (9.71 
4.91 (5.0£ 
6.63 (6.7^ 
4.38 (4.5: 



0.06 
9.49 
5.05 
6.69 
4.49 



0.18 
9.11 
5.14 
6.52 
4.35 



0.03 
9.33 
5.09 
6.52 
4.50 



Meschel et aI.El 
Umemoto et al.EH 



ever, the present MEAM-and other interatomic poten- 
tials in the literature-could not reproduce the low value 
of C44 reported by DFT.ESl 



TABLE VII. Single-crystal elastic moduli, Cxy (GPa) of ce- 
mentite. The predictions from the present work (MEAM) 
are comp ared with those from DFT^ and other interatomic 
potentialpEH. 



Property 


DFT 


MEAM 


Ruda 


Henriksson 
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134 
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135 
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123 


134 
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B. Polycrystalline elastic properties 

Theoretical upper and lower bounds for the polycrys- 
talline bulk modulus (B) and shear modulus (G) were 
calculated using the single-crystal elasti c con stants ac- 
cording to methods by Reuss and Voigtl^SEU xhe true 
polycry stallin e B and G were then estimated using Hill's 
average 1^2'^ Young's modulus (F) and Poisson's ratio 
(v) were calculated by usingPS 

Y ^9BG/{3B + G) (1) 
iy={W/2-G)/{3B + G) (2) 

Polycrystalline elastic moduli predicted by our MEAM 
potential are presented in Table |VIII[ in comparison with 
DFT, experiment, and interatomic potentials by Ruda et 
alf^ and Henriksson et al!^ The elastic constants pre- 
dicted by DFT are in excellent agreement with exper- 
iment. The present MEAM gives the best agreement 
with experiment among the three interatomic potentials 
for B; however, the present MEAM prediction of i' is 
significantly lower than its experimental value. Ruda's 
potential predicts the best agreement with experiment 
for G and v. 



TABLE VIII. Polycrystalline cementite properties: bulk 
modulus B (GPa), shear modulus G (GPa), Young's modu- 
lus Y (GPa), and Poisson's ratio u. The predictions from the 
present work (MEAM) are comp ared w ith those from DFT^^, 
and other interatomic potential^^^^. Experimental values 
are given in parentheses. 

Property DFT (ExpU MEAM Ruda Henriksson 
B 224 (175±^ 174 183 234 

G 72 (7^^ 146 69 114 

Y 194 (177[^ 19^20(Q 343 184 293 
ly 0.36 (0.3^ 0.17 0.33 0.29 

Scott et aI.Eil 
^ Laszio et al.ESl 
" Mizubayashi et alES 

Umemoto et all^ 



C. Surface energies 

Calculations were performed on (001), (010), and (100) 
surfaces to determine the surface formation energy. Ta- 
ble |IX| compares the surface formation energies of the 
present MEAM to DFT and the interatomic potential by 
Ruda et al.^ The atoms near the surfaces are fully re- 
laxed to allow reconstruction if necessary. The predicted 
surface energies have the same order of magnitude as 
DFT results. However, the present MEAM gives a wrong 
order of stability among the three surfaces. 



TABLE IX. Surface energies (J/m^) of cementite. The pre- 
dictions from the present work (MEAM) are compared with 
those from DFT^^ and other interatomic potentiali'^. 



Surface 


DFT 


MEAM 


Ruda 


(001) 


2.05 


2.30 


1.96 


(010) 


2.26 


1.81 


2.00 


(100) 


2.47 


1.79 


2.34 



V. THERMAL PROPERTIES OF CEMENTITE 

A. Thermal stability of cementite 

The stability of cementite at high temperatures was 
investigated through molecular dynamics (MD) simula- 
tions in a canonical (NVT) ensemble at 300 K and 800 K. 
At the end of these MD simulations, cementite retained 
its crystalline structure, affirming its stability at high 
temperatures. The present Fe-C MEAM potential was 
also used to predict several thermal properties of cemen- 
tite. In this section, we present calculations for predict- 
ing melting temperature and variation of specific heat 
and volume of cementite with respect to temperature. 



B. Melting temperature simulation 

The melting temperature of cementite is not well de- 
fined due to its instability at high temperatures. The Fe- 
C phase diagram indicates a eutectic point at 1420 
where liquid consisting of Fe and C solidifies to form 
austenite and cementite crystals. For the purpose of this 
calculation, we considered the melting temperature of ce- 
mentite to be the temperature when cementite loses its 
crystal structure and becomes a random collection of Fe 
and C atoms similar to the eutectic point. The melting 
temperature calculation can be done using a single-phase 
simulation box. However, the single phase method gen- 
erally overestimates the melting temperature due to the 
lack of the interface effects. I^Sl To avoid this superheating 
problem and predict the melting temperature more accu- 
rately, we used a two-phase simulation box that contains 
both solid and liquid phases. 



1. Preparation of two-phase simulation box 

We performed two-phase simulations (TPS) in the 
isothermal-isobaric (NPT) ensemble to determine the 
melting temperature of cementite. The simulation box 
contained both solid and liquid phases of cementite. First 
a supercell containing 14 x 7 x 7 unit cells of cementite 
(10976 atoms) was heated via MD runs in the NPT en- 
semble with T = 1200 K and P = 0. Next, one half of 
the atoms in the supercell were fixed in their positions 
and MD runs were carried out for the other half in the 
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FIG. 8. (Color online) Snapshots of the two-phase MD simulation in the NPT ensemble with T = 1420 K and 
P = Q. Red spheres are Fe atoms and blue spheres are C atoms, (a) Initial state of the simulation box, which contains 
both liquid and solid phases of cementite. (b) Intermediate state of the simulation box at 16 ns, as the solid phase prop- 
agates to the liquid phase, (c) Final state of the simulation box at 32 ns, when the entire system has turned into the solid phase. 




FIG. 9. (Color online) Snapshots of the two-phase MD simulation in the NPT ensemble with T = 1430 K and P = 0. Red 
spheres are Fe atoms and blue spheres are C atoms, (a) Initial state of the simulation box, which contains both liquid and solid 
phases of cementite. (b) Intermediate state of the simulation box at 20 ns, as the liquid phase propagates to the solid phase, 
(c) Final state of the simulation box at 30 ns, when the entire system has turned into the liquid phase. 



NPT ensemble with a sufficiently high temperature (such 
as T = 4000 K) and P = to create a liquid phase. The 
resulting supercell was then subjected to MD runs in the 
NPT ensemble with T = 1500 K (close to T 1200 K 
but still higher than the expected melting temperature) 
and P = 0, still keeping the same half of the atoms fixed. 
The result of this process was a supercell containing solid 
cementite at 1200 K in one half, and liquid cementite at 
1500 K in the other. This ensures a minimum difference 
of stress between atoms in liquid and solid phases of the 
supercell. This supercell was then used in the simulations 
of solidification and melting of cementite. 



2. Two-phase simulation 

The two-phase supercell prepared in the previous sec- 
tion was heated by MD runs in the NPT ensemble where 
the temperature T was increased from 1000 K to 1600 K 
in 100 K intervals. Each system ran for 1.6 ns of simula- 
tion time. The phase change of the two-phase simulation 
box was visually monitored. At 1400 K, the solid phase 
of the simulation box progressed to occupy the entire 
box, and at 1500 K the liquid phase of the simulation 
box progressed to occupy the entire box. Next, the ini- 
tial two-phase simulation box was heated from 1400 K to 



1500 K in 10 K intervals using NPT MD runs. Each sys- 
tem was equilibrated for at least 5x10^ time steps, where 
each time step was 2 fs, totaling to 10 ns. The final state 
of the system was visually inspected. If the final state 
appeared to have both liquid and solid phases, more MD 
runs were performed until the final state of the super- 
cell contained only one phase. Some systems required as 
much as 32 ns of MD runs to arrive at a single phase. 
The transformation of the two-phase simulation box to 
a one-phase simulation box near the predicted melting 
temperature is presented in Fig. [8] and Fig. [9] The total 
energy, volume, and pressure of the systems were deter- 
mined through averaging the values of the final 40, 000 
time steps (80 ps) of each simulation. 

In Fig. [To] we plot the total energy, volume, specific 
heat, and the derivative of volume as functions of tem- 
perature. Experimental data for specific heat and vol- 
ume are not available for the 1400-1500 K temperature 
range. Available experimental data are the heat capacity 
of 3.6 fcs/atom at 1023 K,^^ which converts to 3.4 eV/K 
for the current simulation, and the experimental vol- 
ume of 10 A^/atom at 1070 K.S51 Specific heat and vol- 
ume determined by Dick et al. from the first-principles 
calculation^^ done on the solid phase of cementite are 
included for comparison in Figs. 10 b) and (c). Since 



Dick and coworkers used a single-phase simulation box. 
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their simulation clearly shows superheating causing the 
melting temperature to be overestimated. This can be 
attributed to the absence of the solid-liquid interface in 
single phase simulations. In Fig. [lOjjc) the specific heat 
shows a peak between 1420 K and 1430 K. Therefore we 
assign 1425 ± 5 K as the melting temperature of cemen- 
tite. This is in good agreement with the experimental 
eutectic point of 1420 K.— 

VI. SUMMARY AND CONCLUSION 

We investigated the properties of cementite using an 
interatomic potential developed within the MEAM for- 
malism. Previously developed single element interatomic 
potentials for Fe and C were used to develop the Fe-C al- 
loy MEAM potential. The single element potential for C 
predicts graphite and diamond as minimum energy struc- 
tures with almost degenerate energies. MEAM poten- 
tials for pure elements predict the heat of formation, bulk 
moduli, and lattice constants of Fe and C in their natural 
crystal structures in good agreement with experimental 
data. The alloy potential for the Fe-C system was devel- 
oped to reproduce heat of formation of Fe-C alloys in Bi, 
Li2, and cementite crystal structures, as well as the inter- 
stitial energies of C in bcc Fe. The Fe-C potential was val- 
idated by investigating the energy variation with respect 
to volume of Fe-C alloys in Bi and L12 structures. The 
validated potential was used to predict structural, elastic. 



and thermal properties of cementite. Structural proper- 
ties tested included the heat of formation, the equilibrium 
lattice constants, the equilibrium volume, and the energy 
variation with respect to volume. MEAM predictions are 
in good agreement with DFT and experiment. The nine 
single crystal elastic constants were calculated and used 
to estimate polycrystalline bulk modulus, shear modu- 
lus. Young's modulus, and Poisson's ratio of cementite. 
Surface energies for (001), (010), and (100) surfaces were 
also calculated and compared. Melting temperature and 
the variation of specific heat and volume with respect to 
temperature were predicted by two-phase (solid/liquid) 
MD simulations. The present MEAM potential predicted 
the melting temperature of cementite to be 1425 ± 5 K, 
which is in good agreement with the eutectic point of 
1420 K given in the Fe-C phase diagram. 
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